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ABSTRACT 

The  long-range  and  short-range  order  parameters  are  computed 
for  the  Ising  lattice  using  a  Monte  Carlo  sampling  scheme.  The 
square  lattice,  the  simple  cubic  lattice  and  the  body-centered  cubic 
lattice  are  considered.   In  the  three  dimensional  calculations  both 
the  antiferromagnetic  and  ferromagnetic  cases  are  considered  as  well 
as  the  coupling  to  an  external  magnetic  field  of  various  strengths. 
Good  agreement  is  found  where  the  results  can  be  compared  with  the 
exact  two  dimensional  results  and  in  the  three  dimensional  case  the 
results  agree  well  with  those  obtained  from  series  approximations  in 
the  regions  where  the  series  approximations  are  valid.  The  present 
method  appears  to  give  good  results  for  the  short-range  order  even 
very  close  to  the  critical  temperature,  but  in  this  neighborhood  the 
long-range  order  estimate  is  crude.  The  computations  were  performed 
on  the  high-speed  computer  IT  T.I  AC,  located  at  the  University  of 
Illinois . 


I.      INTRODUCTION 

(1) 

In  a  recent   investigation x    '   the  Monte  Carlo  method  was   used  to 

compute  parameters  describing  the  short-range  and  long-range  order  in  a  face- 
centered  cubic  binary  alloy.     That   investigation  was  preceded  by  some  compu- 
tations^   '   of  order  parameters   in  a  two-dimensional  Ising  lattice^    '   as  a  test 
of  the  feasibility  of  the  Monte  Carlo  method  for  this  kind  of  calculation. 
This  early  work  was   quite  successful  and  the  continuation  of  this  work  to  a 
treatment  of  three-dimensional  Ising  lattices   is  the  subject  of  the  present 
paper.     Since  the  early  work  on  the  two-dimensional  lattice  has  not  been 
previously  reported  in  detail,    it  has  been  included  in  the  present  discussion. 
Two  three-dimensional  lattices   are  treated  in  the  present  works      the  simple 
cubic  and  the  body-centered  cubic.     This  treatment  includes  both  ferromagnetic 
and  antiferromagnetic  coupling,   and  coupling  to  an  external  magnetic  field. 
Parameters  describing  the  short-range  order  and  the   long-range  order  have  been 

computed . 

(k) 

The  method  used  here  was  first  used  by  Metropolis  and  others v  '   to 

treat  the  two-dimensional  hard  sphere  gas,  and  it  has  been  used  subsequently 

by  others  ^     '   for  further  computations  on  the  equation  of  state  of  gases. 

The  essence  of  this  method  can  be  described  briefly  as  follows.  A  sequence  of 

configuration  states  for  the  system  is  developed  using  transition  probabilities 

p. .,  where  p. .  gives  the  probability  that  state  i  will  be  followed  immediately 
J- j        J-  j 

by  state  j.  These  transition  probabilities  are  chosen  to  make  the  distribution 
of  states  in  the  sequence  tend  toward  a  Boltzmann  distribution  as  the  number  of 
states  in  the  sequence  increases.  At  some  point  the  sequence  is  truncated  and, 
neglecting  some  of  the  initial  states  in  the  chain,  the  states  of  the  truncated 
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sequence  are  used  as  an  ensemble  to  estimate  the  average  value  of  certain  system 
parameters;  in  the  present  case,  estimates  of  the  average  value  of  the  order  para- 
meters are  computed.  It  seems  appropriate  to  describe  this  approach  as  a 
"mathematical  experiment"  because  it  is  somewhat  analogous  to  observing  the 
parameters  directly  in  the  real  physical  system,  as  in  a  physical  experiment. 
In  the  latter  case  nature  provides  the  averaging,  whereas  in  the  mathematical 
experiment  this  is  simulated  on  a  model.  It  should  be  quickly  pointed  out, 
however,  that  the  kinetics  associated  with  the  mathematical  experiment  do  not 
necessarily  represent  those  of  the  real  system;  they  may  represent  the  real 
system  kinetics  to  some  degree,  but  it  is  not  essential  to  the  method.  This 
approach  can  provide  a  very  good  physical  picture  of  the  microscopic  character 
of  the  system  and  it  is  therefore  capable  of  providing  new  insights  to  the 
problem  which  might  be  very  difficult  to  obtain  from  a  more  conventional, 
analytical  approach  in  which  the  system  is  represented  in  a  comparatively 
abstract  fashion. 

Without  a  high-speed  computing  machine  this  approach  would  not  be 
feasible.  It  is  not  surprising  therefore  that  interest  in  this  method  has 
increased  as  these  machines  have  become  more  available,  and  it  is  to  be  expected 
that  this  interest  will  continue  to  grow  as  the  capabilities  of  these  machines 
grow.  The  ILIIAC,  a  high-speed  computing  machine  located  at  the  Digital  Computer 
Laboratory  of  the  University  of  Illinois,  was  used  to  perform  the  computations 
described  in  the  present  work. 
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II.   GENERAL  REMARKS  ON  THE  APPLICATION  OF  THE  MONTE  CARLO 
METHOD  TO  THE  ISING  LATTICE  PROBLEM 

The  fundamental  ideas  of  the  method  which  was  used  in  these  computations 
have  "been  discussed  in  the  references  cited  in  the  introduction,  especially 
references  (l),  (k)   and  (6).  A  familiarity  with  these  ideas  will  be  assumed 
and  here  our  attention  will  be  directed  at  their  application  to  the  Ising 
lattice  problem. 

Each  site  of  the  Ising  lattice  has  an  associated,  two-valued,  spin 
coordinate  P-(k)  for  the  k —  site$  n(k)  =  +  1  or  -  1,  corresponding  to  the  two 
allowed  orientations  of  the  spin  on  the  k —  site.  The  i —  configuration  state 
of  a  lattice  of  N  sites  is  completely  descrihed  hy  the  N- component  vector  a  , 
whose  components  are  the  N  spin  coordinates,  |i.(k).  The  energy  of  this  spin 
array  is  assumed  to  he  due  to  nearest-neighbor  interactions  and  the  interaction 
with  an  external  magnetic  field.  In  particular,  the  energy  of  the  i —  configu- 
ration is  given  by 

E  =  -J  E(l)  H.(k)  u.(k')  +  Hi  n  (k)  (1) 

1     k^k     X     X         k  X     > 

where  the  first  sum  extends  over  all  pairs  of  sites  in  the  lattice  such  that 
k  and  k'  are  nearest-neighbors  and  the  second  sum  extends  over  all  sites  in 
the  lattice.  It  will  be  recognized  that  the  change  in  the  coupling  energy  for 
a  nearest-neighbor  pair  going  from  a  state  of  parallelism  to  a  state  of  anti- 
parallelism  is  2J:  i.e.,  l  (tl)  -  C   (tf)  =  2J.  The  array  is  ferromagnetic 
if  J  is  positive  and  antiferromagnetic  if  J  is  negative.  Referring  again  to 
Eq.  (l),  it  will  be  observed  that  the  change  in  the  energy  arising  from  the 
external  field  coupling  is  2H  when  a  spin  changes  from  a  state  of  parallelism 
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with  the  field  to  a  state  of  ant i parallelism  with  the  field;  i0e0, 
C   (t  )  "  C  (i)  -   2H.   Since  it  is  convenient  and  customary  to  use  para- 
meters which  are  the  ratios  of  the  coupling  energies  to  kT,  where  k  is  the 
Boltzmann  constant  and  T  is  the  absolute  temperature,  we  define  K  =  j/kT  and 
L  =  H/kT. 

The  order  in  the  system,,  at  equilibrium,  is  computed  as  a  function 
of  K  and  L.  To  describe  the  short-range  order,  the  parameters  f.(j)  are  used 

where  f . ( j)  is  the  fraction  of  j — -neighbor  sites  which  are  occupied  by  an 

th 
antiparallel  pair  of  spins  in  the  i—  configurations 


f*(J)  =  A  |f 


1  -  u.(k)  u,(k') 


(2) 


th 
where  cc(j)  is  the  number  of  j — -neighbors  of  a  lattice  site,  and  the  summation 

th 
extends  over  all  pairs  of  sites  k  and  k'  where  k  and  k'  are  j — -neighbors. 

th 
Using  Boltzmann  statistics,  the  average  value  of  this  parameter  for  j — - 

neighbors  is  given  by  the  usual  formula 


-E.  /,„ 

tin  =  s  f,(j)  -~—  .  (3) 

i 
where  the  summations  extend  over  all  configuration  states  of  the  system.   In 

the  Monte  Carlo  estimation  of  this  quantity  the  expression  on  the  right  is 

replaced  by  the  average  over  a  small  sample  of  f.(j)'s  drawn  from  a  distri- 

th 
bution  in  which  the  i —  configuration  state  tends  to  be  proportional  to 

e   x'   o  This  proportionality  holds  strictly  only  in  the  limit  as  the 

sequential  process  of  generating  new  configurations  is  continued  indefinitely. 

However,  the  proportionality  usually  holds  with  sufficient  accuracy  for  worth- 


while  results  when  the   sequence   is  truncated  in  order  to  keep  the   computing 
time  within  reasonable  hounds.     For  the  two-dimensional  square   lattice  and 
the  body- centered  cubic   lattice  the  average  of  the  first-neighbor  order  para- 
meter,  f(l),    is   estimated.      In  the  simple- cubic   lattice  both  f(l)    and  f(2) 
are  estimated. 

The   long-range  order  of  the   i — -  configuration  state   is   described  by 
the  parameter  S.,  where 


•i-i 


E  |i,  (k) 


k    x 


(M 


the  summation  extending  over  all  sites.  The  average  value  of  S.,  denoted  by 
S,  is  given  by  the  formula  analogous  to  Eq.  (3)  and  the  Monte  Carlo  estimate 
of  S  is  likewise  obtained  from  sampling  as  indicated  above  for  f(j).  Estimates 
of  S  have  been  obtained  for  the  two  three-dimensional  lattices „      It  is  known 
that  for  an  infinite  lattice  and  L  =  0,  S  vanishes  at  a  critical  value  of  K, 
denoted  by  K  ,  and  that  for  K  <  K  ,  S  remains  equal  to  zero.  The  finite  system 
used  in  these  computations  can  only  be  expected  to  approximate  this  behavior, 
and  one  expects  a  rapid  decrease  in  S  in  the  neighborhood  of  K  ,  but  S  will 

remain  non-zero,  though  small,  even  when  K<<  K  .  When  the  external  magnetic 

(9) 
field  is  zero,  S  can  be  identified  as  the  spontaneous  magnetization  per  spin. 

For  a  lattice  composed  of  an  infinite  number  of  spins  a  very  small  positive 
magnetic  field  removes  all  states  of  the  lattice  from  the  ensemble  average  for 
which  the  total  spin  is  negative,  but  with  a  finite  lattice  this  is  not  strictly 
true.  This  approximation  to  the  spontaneous  magnetization  of  an  infinite  array 
will  be  worst  near  the  critical  temperature  where  the  difference  in  behavior  be- 
tween the  infinite  system  and  the  finite  system  becomes  particularly  important 0 


The  procedure  for  generating  the  ensemble  of  configurations  over  which 
the  averaging  is  to  he  performed  is  very  similar  to  the  one  described  in  refer- 
ence (l).  We  are  given  a  lattice  of  N  sites  which  is  in  a  particular  configura- 
tion state,  say  i;  hence  M-.(k)  is  known  for  all  k.  A  site  of  the  lattice  is 
selected  and  the  change  in  energy,  7\Ej  which  would  accompany  a  reversal  in 
orientation  of  the  spin  on  that  site  is  computed.   Tf  S\K  <  0,  then  the  spin 

is  reversed  and  if  /ATH  >  0,  then  the  spin  orientation  is  reversed  with  proba- 

-^\E/kT 
bility  e    '   .  The  latter  process  is  done  by  generating  a  pseudo-random 

number,  $  ,  from  a  uniform  distribution  on  the  interval  (0,  l)  and  if 

<r^    -^JS/kT 
_5  <  e    '   ,  the  spin  orientation  is  reversed;  otherwise  it  is  not  reversed. 

Next,  another  site  is  selected  and  the  above  process  is  repeated.  The  sites 
were  numbered  and  selected  sequentially  and  when  all  sites  of  the  lattice  had 
been  treated  once,  as  described  above,  an  iteration  of  the  calculation  was 
said  to  be  completed.  At  the  end  of  an  iteration  the  values  of  the  order  para- 
meters f.(j)  and  S.  were  recorded;  the  details  of  this  step  varied  somewhat  in 
the  different  calculations  and  they  will  be  described  more  fully  in  the  following 
sections.  The  iterations  continued  until  a  condition  for  stopping  was  satisfied 
at  which  point  the  last  R  samples  which  were  generated  were  used  to  compute  the 
averages: 

f(i)  =  "|  Z  f.(i), 

(i)  x 


f(2)  -  -|   I  f.(2),  (5) 

R  (i)   X 

R  (i)   * 
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where  in  each  case  the  sum  extends  over  the  last  R  samples.  In  the  three 
dimensional  computations  the  stopping  condition  was  based  on  a  comparison  of 
results  from  two  independent  computations,  similar  to  that  described  in 
reference  (l);  the  details  will  be  presented  later. 

It  is  not  difficult  to  see  that  this  process  satisfies  the  necessary 

(1) 

and  sufficient  conditions^   for  producing  a  sequence  of  configurations  which 

tends  toward  a  Boltzmann  distribution  as  the  length  of  the  sequence  tends 
toward  infinity.  The  ergodic  condition  is  satisfied"  Any  pair  of  configu- 
rations can  be  linked  by  at  least  one  finite  sequence  of  configurations  in 

which  one  spin  at  a  time  is  reversed,  and  the  probability  for  this  reversal 

N 
is  non-zero.  The  other  condition  is  also  satisfied?  The  2  component  proba- 
bility vector  j  with  components 

* 1  "y  'Ei/KT  (1  =  *>    2>     ■■•>    2)  (6) 

i 

is  an  eigenvector,  with  eigenvalue  unity,  of  the  stochastic  matrix  P,  whose 

N 
elements,  p  .,  are  the  transition  probabilities  linking  the  2  configuration 

states  of  the  lattice  in  a  complete  iteration.  This  can  be  seen  as  follows. 

The  stochastic  matrix  P  may  be  regarded  as  the  product  of  N  matrices 

P(l),  P(2)„.„P(N)  where  P(k)  has  components  p.  .(k)  which  are  the  probabilities 

for  a  transition  from  state  i  to  j  by  a  reversal  in  the  orientation  of  the  spin 

on  site  k.  For  given  i  there  is  exactly  one  j,  say  j',  for  which  p„  .(k)  =/  0 

and  p„.(k)  =  1  -  p.  .,(k).  If  E.  >  E.,,  we  have  p.  „ ,  (k)  =  1  and  if  E.  <  E.f,  we 

.=  ("p1       =.1?  ^ 

have  p,  .  ,(k)  =  e  ^  j'   i'/KT°   Similarly,  for  given  j  there  is  exactly  one  i, 
say  i'  for  which  p.  .  (k)  ^   0,  where  i"  ^   j.   In  short,  there  will  be  just  one 


or  two  non-zero  elements  in  every  row  and  column  of  P(k) .     Consider  the 
product 


fp(k)  =$ 


and  in  particular 


i     i        J  d 


(7) 


Let  p.,  .(k)   be  the  non-zero  off-diagonal  element  in  the   j —  column  of  P(k) 

1   J  ■(BI'~E1)/kT 

If  E.  ,  >  E.,   then  p.  ,  .  =   1  and  p  .  .  =   1  -e  J  '        and  it  follows  from 

Eqs.    (6)    and   (7)  that 


0.  =  Ae 


/kT 


-(E..-E.) 


1-e 


dVKP 


Ae 


'Ej/kT 


hence 


where 


4.  =  Ae"^/^ 


A-i  a  ^i/ia. 


If  E.  ,  <  E.,   then  p.  ,  . 
Eqs„    (6)    and   (7)  that 


=  e 


(VEi?)/kT 


and  p  .  .  =  0  and  it  follows   from 


J 


(E-E,,)/kT 


J      i 


-E. 


Ae  , 


/kT 


hence 


Qi.  =  Ae 


■E 


j/kT 
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It  follows  that  $  =  'Jr  ,  hence  y   is  an  engenvector  of  P(k)  with  eigenvalue 

unity  and  since  this  is  true  for  all  k  it  is  true  for  the  product 

P  =  P(l)  P(2)"*P(N).   It  is  to  be  noted  that  this  result  is  independent  of 

the  order  in  which  the  sites  are  numbered.   This  result  is  also  true  for 

1  N 
Q  =  —  nZ-.  P(k),  which  is  the  transition  matrix  when  a  site  is  selected  at 
N  k=l     ' 

random  and  the  same  process  performed  on  it.   The  details  of  the  individual 
computations  follow. 

III.   THE  SQUARE  LATTICE 

The  two  critical  items  affecting  the  practicality  of  the  method, 
namely  the  rate  of  convergence  of  the  generated  ensemble  to  a  Boltzmann 
distribution  and  the  size  of  the  lattice  needed  for  worthwhile  results,  were 
investigated  using  this  model.  The  short-range  order  parameter  f (l)  was 
computed  and  compared  with  the  exact  value  for  an  infinite  system. *   '  The 

_OV" 

computations  were  made  at  different  values  of  the  parameter  x  =  e    and  with 
zero  external  field,  L  =  0.   Periodic  boundary  conditions  were  imposed  by 
linking  the  left  edge  of  the  lattice  to  the  right  edge  by  nearest- neighbor 
bonds  and  similarly  the  top  edge  was  linked  to  the  bottom  edge. 

The  lattice  configuration  was  represented  by  binary  numbers  in 
the  computer,  with  each  binary  digit  corresponding  to  a  lattice  site  and 
the  value  of  that  digit  representing  the  orientation  of  the  spin  on  the 
site.   In  each  iteration  of  the  computation  the  sites  were  selected  for 
consideration  systematically.  The  successive  sites  along  a  row  were  con- 
sidered until  the  end  of  the  row  was  reached  and  then  the  sites  in  the  next 
adjacent  row  were  considered  until  finally  all  N  sites  had  been  considered. 
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At  the  completion  of  one  iteration  the  value  of  f . (l)  was  recorded. 
To  examine  the  hehavior  of  estimates  of  f (l)  as  it  is  computed  at  different 
points  along  the  chain,  the  sequence  of  values  of  f.(l)  was  broken  into  groups 
of  128  and  the  average  and  standard  deviation  computed  for  each  of  the  groups  1 
the  sequence  of  averages  thus  generated  will  be  denoted  by  f(l,  1),  f (l,  2), 

f(l,  3)---. 

The  computations  were  performed  for  three  different  lattice  sizes  1 
10  sites  on  an  edge,  20  sites  on  an  edge  and  37  sites  on  an  edge.  The  com- 
puting time  to  complete  one  iteration  was  approximately  2  seconds  for  the 
10  x  10  lattice  and  this  time  varies  linearly  with  the  number  of  sites  in  the 
lattice. 

In  order  to  observe  the  effect  of  the  choice  of  the  initial  configu- 
ration on  these  results,  all  of  the  computations  were  performed  twice,  using 
quite  different  starting  conditions  for  the  two  cases.  In  one  case  the 
initial  configuration  was  one  of  complete  order  with  all  spins  up|  that  is, 
n(k)  =  1  for  all  k.  In  the  other  case  the  initial  configuration  was  highly 
disordered^  this  configuration  was  generated  by  assigning  the  values  of  M-(k) 
such  that  there  would  be  equal  probability  for  up  and  down  spins. 

The  results  of  the  calculation  of  f (l9   j)  are  shown  in  Table  I. 
From  the  exact  treatment  of  the  square  Ising  lattice  it  is  known  that  there 
is  a  second  order  phase  transition,  located  at  x  =  0.^1^-2,  at  which  the  con- 
figurational  specific  heat  becomes  logarithmically  infinite.  Since  the 
specific  heat  is  proportional  to  the  variance  of  the  energy 


d  <E> 


oC<[E  -<E>]2>        , 


d  T 
it  is  not  surprising  that  the  strongest  fluctuations  and  largest  errors  in 
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the  results  occur  in  the  neighborhood  of  x  =  0o^3»  It  is  interesting  to  note 
that  standard  deviations  vary  proportionally  to  l/-\/N~as  is  to  be  expected  for 
purely  statistical  reasons „   It  is  satisfying  to  find  that  a  prohibitively 
large  lattice,  and  the  37  x  37  lattice  approaches  this,  is  not  really  needed 
for  a  good  estimate  of  f(l)»  In  fact  the  difference  in  accuracy  between  the 
two  larger  systems  can  hardly  be  called  significant „   Hence,  it  appears  that 
beyond  the  20  x  20  lattice  a  large  increase  in  N,  and  consequently  a  large 
increase  in  computing  time,  would  be  needed  for  a  small  increase  in  accuracy. 
It  also  appears  that  the  effects  of  the  initial  configuration  are  lost  rather 
quickly  and  in  fact  the  average  taken  over  the  first  sample  of  128  configurations 
is  frequently  in  good  agreement  with  the  exact  values  this  is  particularly  true 
for  the  system  which  started  from  a  completely  ordered  configuration  It  should 
be  pointed  out  here  that  the  very  first  configuration,  the  completely  ordered 
one  or  the  disordered  one,  is  not  included  in  the  averaging. 


IV.  THE  SIMPLE  CUBIC  LATTICE 

The  ILLIAC  program  which  was  prepared  to  do  the  computations  on  the 
simple  cubic  lattice  is  more  elaborate  than  the  one  which  was  prepared  for  the 
square  lattice  computations.  An  important  part  of  the  present  program  is  a 
test  to  determine  a  point  in  the  sequence  of  configurations  at  which  the 
computation  of  the  ensemble  averages  is  to  commence .  This  test,  called  the 
convergence  test,  resembles  a  similar  test  used  in  the  work  described  in 
reference  (l),  and  its  present  application  is  explained  below. 
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The  problem  of  obtaining  a  reasonably  accurate  result  without  using 
enormous  amounts  of  computing  time  depends  partly,  as  mentioned  earlier,  on 
the  rate  of  convergence  of  the  generated  ensemble  to  a  Boltzmann  distribution. 
The  work  on  the  square  lattice  has  shown  that  in  many  cases  the  convergence 
is  quite  rapid.  With  the  square  lattice  the  exact  solution  provides  a  guide 
for  checking  the  ensemble  averages,  but  when  the  exact  solution  is  not  known 
a  rule  must  be  made  for  selecting  the  point  at  which  the  averaging  may  commence, 
The  problem  of  constructing  such  a  rule  is  tricky.  One  might  consider  the 
successive  values  of  a  parameter,  such  as  the  short-range  order  parameter,  and 
commence  averaging  where  this  sequence  appears  to  be  steady  in  some  sense. 
This  is  dangerous  for  there  may  be  two,  or  more,  sets  of  states  each  of  which 
has  a  relatively  high  probability  of  occuring,  but  which  are  linked  by  small 
transition  probabilities.   Such  a  situation  would  result  in  fairly  steady 
sequential  values  for  certain  parameters  except  that  here  and  there  the 
apparent  equilibrium  value  of  the  parameter  might  change  abruptly,,   It  may 
happen  that  these  abrupt  changes  occur  so  infrequently  that  they  would  not 
even  occur  in  a  very  long  (on  the  computational  time  scale)  computation. 
In  this  instance  an  average  taken  over  the  apparently  steady  sequence  might 
give  a  very  incorrect  result 0     Since  the  region  in  phase  space  which  con- 
tributes significantly  to  the  ensemble  average  becomes  broader  near  a  phase 
transition  it  is  to  be  expected  that  this  phenomenon  is  likely  to  occur  in 
such  a  region,,  It  is  true  that  no  such  difficulties  were  ever  apparent  in 
the  two  dimensional  lattice  computations  but  similar  difficulties  have  been 
encountered  in  the  hard  sphere  equation  of  state  computations. 

The  rule  which  has  been  adopted  for  the  present  computations  is 
characterized  by  the  fact  that  two  statistically  independent  sets  of  results 
are  developed  and  the  point  of  convergence  is  established  when  these  results 
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agree  within  a  given  margin  of  error.  Two  distinct  lattices  are  used  to  develop 
two  statistically  independent  sequences  of  configurations.   One  of  the  sequences 
starts  from  a  configuration  of  complete  order,  while  the  other  starts  from  a 
disordered  configuration,  just  like  the  two  initial  configurations  discussed  in 
the  square  lattice  computations.  The  difference  between  the  present  case  and 
the  former  is  that  now  the  two  sequences  are  developed  simultaneously  and  hence 
may  he  compared,  one  with  the  other,  as  the  two  sequences  are  developed.  We 
denote  the  sequence  starting  from  a  configuration  of  complete  order  as  the  low 
temperature  (IT)  sequence.   Correspondingly,  the  sequence  starting  from  a  dis- 
ordered configuration  is  called  the  high  temperature  (HT)  sequence.  The  gener- 
ation of  the  ensemble  is  divided  into  three  stages.   In  the  first  stage  a  LT 
sequence  of  NL  configurations  and  a  HT  sequence  of  Mn  configurations  is  devel- 
oped. Three  parameters  associated  with  each  configuration  are  held  in  the 
computer  stores   f.(l),  f.(2)  and  S.„   In  the  second  stage  the  average  of  f . (l) 
taken  over  the  last  JYL  configurations  is  computed  for  each  sequence  and  similarly 
for  f.(2);  these  are  designated  p(l)l  JJJ;9      \t  (l)J  m,      ^(2)1  ^  and  1^(2)1^. 
The  convergence  test  is  passed  when  the  following  two  inequalities  are  satisfied 
for  the  first  time; 


<*i 


(8a) 


[f(2)]M  +  pjjm 


<£. 
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where  t,  and  £   are  small  positive  numbers.   If  both  conditions  are  not 
satisfied  a  new  configuration  is  generated  for  the  IT  sequence,  and  for  the 
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HT  sequence.  The  oldest  information  in  the  sequence,  namely  that  pertaining  to 
the  configuration  which  occurred  M_  iterations  before  the  present  point,  is 
thrown  out  to  make  space  for  the  information  on  the  new  configuration.  Thus, 
the  order  parameters  for  the  last  M  configurations  in  each  chain  are  retained., 
As  each  new  configuration  is  added  to  the  UP   sequence  and  to  the  HT  sequence 
the  convergence  test  is  repeated.  When  the  convergence  test  is  finally  passed 
the  third  stage  of  the  computations  begins.  All  of  the  values  for  f . (l)  in 
the  two  sequences  are  collected  into  one  sum,  E  f.(l),  and  similarly  for  f.(2), 
and  S.o  As  each  new  configuration  is  generated  in  each  sequence  the  values  for 
f.(l),  f.(2)  and  S.  are  added  to  the  corresponding  sum.  When  ^\M  new  configu- 
rations in  each  sequence  have  been  generated  then  the  sums  are  divided  by 
2(Mn  +Z^M)  =  R,  the  number  of  configurations  in  the  ensemble,  to  obtain  the 
final  estimates  of  the  quantities  f(l),  f(2)  and  S„   In  the  third  stage  of  the 
computations  the  sum  of  squares  of  each  of  the  parameters  is  also  generated  in 
order  to  calculate  the  standard  deviations. 

This  process  does  not  insure  against  the  possibility  of  obtaining 
an  erroneous  answer  because  of  the  chains  being  trapped  in  a  set  of  meta- 
stable  states.  It  can.be  expected,  however,  that  the  chance  of  detecting  such 
a  situation  is  better  than  it  would  be  if  only  one  sequence  was  considered. 
Of  course,  if  memory  space  and  computing  time  permit,  then  one  can  extend  this 
method  to  include  a  still  larger  number  of  independent  chains. 

The  model  for  almost  all  of  the  computations  had  8  sites  on  an  edge, 
and  thus  contained  a  total  of  512  sites.   Some  computations  were  done  on  a 
model  with  16  sites  on  an  edge  but  because  of  the  large  amounts  of  computing 
time  required  this  work  was  rather  limited.  The  amount  of  time  required  to 
complete  one  iteration  for  the  8x8x8  array  was  about  6  seconds  and  it  was 
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approximately  eight  times  this  for  the  l6  x  l6  x  l6  array.  As  with  the  square 
lattice,  periodic  boundary  conditions  were  always  imposed.  The  sites  were 
selected  systematically  for  detailed  consideration  in  an  analogous  fashion  to 
the  method  used  with  the  square  lattice:   successive  sites  in  a  row  were 
treated,  then  successive  rows  and  finally  successive  planes. 

The  results  of  this  computation  are  compiled  in  Table  II.   Order 
parameters  S,  f(l)  and  f(2)  are  shown  as  functions  of  K  in  Figs.  1,  2  and  3, 
respectively.   In  Table  II  the  numbers  in  the  first  column  are  simply  identi- 
fication numbers  (ID) .  The  next  two  columns  contain  the  energy  parameters  K 
and  L.   In  the  last  three  columns  the  numbers  M_  and  /\M  appearing  there  have 
already  been  defined  and  the  numbers  in  the  column  headed  "total",  are  the 
total  number  of  iterations  performed.   Hence,  in  row  1,  the  values  M_  =  ^>0, 
/\M  =  50,  total  =  156  mean  that  fifty  configurations  were  generated  in  the 
HT  sequence  before  the  convergence  tests  began  (since  M_  =  ^0);   next,  after 
convergence  fifty  more  configurations  were  generated  in  each  sequence  (since 
y\M  =  50)  "to  make  a  total  of  200  configurations  in  the  ensemble;  since  the 
total  number  of  iterations  was  156  it  follows  that  the  convergence  test  was 
passed  upon  completion  of  the  one  hundred  and  sixth  iteration  and  therefore 
the  first  fifty-six  configurations  in  each  sequence  were  discarded.   In  the 
columns  following  K  and  L  the  order  parameters  S,  f (l)  and  f (2)  are  presented. 
The  spread  indicated  for  each  order  parameter  is  the  standard  deviation  of 
the  mean.  The  results  are  listed  in  order  of  decreasing  K  (i.e.  increasing 
temperature  for  fixed  J)  in  groups  in  which  L  is  fixed.   Near  the  end  of  the 
table,  starting  at  ID  =  78,  some  results  are  grouped  together  in  which  L/K  is 
fixed:   note  that  L/K  =  H/j,  the  ratio  of  the  external  field  coupling  energy  to 
the  nearest-neighbor  coupling  energy,  which  is  independent  of  T.   In  the  figures 
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the  results   are  plotted  for  fixed  L  (solid  lines)    and  for  fixed  L/K  (dashed 
lines).     Although  the  order  parameters   can  be  computed  analytically  at  K  =  0, 
the  curves  for  constant  L/K  have  not  "been  extended  through  K  =  0  because  com- 
putations were  not  performed  in  the  region  of  K  =  0  and  an  extrapolation  of 
these  curves   from  the  available  data  did  not   seem  appropriate. 

It  will  be  noted  that  three  different  values  for  M^  appear  in  the 
table:      100,    5°  and  25«      In  ^e  very  first    (chronologically)    computations 
the  large  value  of  M_  was  used  together  with^M  =  300,   but  becuase  of  the 
large  amounts  of  computing  time  being  absorbed  it  was  decided  to   set  Mn  and 
^M  both  at  50.      Still  later,   for  reasons  of  economizing  on  computing  time 
the  still  smaller  values  M-  =  25  and^lM  =  25  were  introduced.      In  the  first 
two  cases  the  convergence  test  parameters  were  given  the  value  0.02: 
£      =     £      =  0.02.      In  the   last  case,    in  an  attempt  to  compensate  for  the 
relatively  small  value  of  NL  =  25  the  parameters  were  given  the  value  0.01: 

(±  =  e>  =  0.01. 

Since  the  size  effect  can  be  expected  to  be  most  significant  in 
the  zero  field  case  in  the  neighborhood  of  the  apparent  critical  temperature, 
some  computations  for  a  16  x  16  x  16  array  were  made  in  this  region:  these 
have  identification  numbers  8,  11  and  13.  When  these  are  compared  against 
the  corresponding  results  for  the  8x8x8  system,  it  will  be  noted  that 
there  is  a  marked  difference  in  the  value  of  S  for  the  two  cases.  The 
differences  in  the  results  obtained  for  the  short-range  order  parameters 
on  the  other  hand  are  relatively  slight.   Hence,  it  appears  that  although 
the  estimate  of  S  is  a  crude  approximation  of  its  value  for  the  infinite 
system  in  this  region,  the  estimates  of  f (l)  and  f (2)  are  rather  good.  In 
the  case  of  L  =  0.125  a  computation  was  made  with  the  large  lattice  (ID  =  33) 
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in  the  region  where  S  can  he  expected  to  he  most  sensitive  to  size  effects. 
Comparison  of  these  results  with  those  for  th  8  x  8  x  8  lattice  shown  that 
the  difference  in  the  results  obtained  for  S,  as  well  as  for  f (l)  and  f (2), 
is  slight o   Hence,  for  this  value  and  higher  values  of  L  the  estimates  of  S 
can  he  expected  to  he  fairly  a  good  approximation  to  the  value  for  an  infinite 
system,, 

In  the  ant i ferromagnetic  region  the  value  L/K  =  -6  is  a  critical  one, 
For  L/K  greater  in  magnitude  than  this  value  the  external  field  coupling  domi- 
nates the  nearest-neighbor  coupling  and  at  low  temperatures  the  system  tends 
to  the  state  in  which  all  spins  are  parallel  to  the  external  field.  For  L/K 
smaller  in  magnitude  than  this  value,  the  nearest-neighbor  coupling  dominates, 
and  at  low  temperatures  the  system  tends  to  the  state  in  which  all  nearest- 
neighbor  spins  are  antiparallelo  The  series  of  computations  at  L/K  =  -k   and 
L/K  sa  -8  illustrate  the  alternate  behavior  in  the  order  parameters  as  the 
parameter  K  tends  to  large  negative  values  (i0e0  as  T  -* 0  for  J  equal  to  a 
negative  constant) „     It  is  interesting  to  notice  that  a  relatively  large 
number  of  iterations  had  to  be  performed  in  the  computation  with  ID  =  92 
before  the  convergence  condition  was  satisfied^  in  computations  85  and 
86  a  similar  situation  is  noticed,,  The  reason  for  this  is  that  we  are  near 
the  critical  magnitude  of  K  while  below  the  critical  magnitude  of  L/k„   In 
the  region  K  <  0,  -6  K  >  L  >  0  there  are  two  states  of  minimum  energy  (the 
two  states  in  which  all  neighbors  are  ant i parallel),  just  as  there  are  on 
the  line  K  >  0,  L  -  0..  Thus,  there  will  again  be  a  critical  value  of  K, 
around  which  configurations  consisting  of  mixtures  of  these  two  states  will 
tend  to  persist o  As  the  line  L/K  -  °6  is  approached,  a  third  configuration 
also  assumes  importance,  that  in  which  all  spins  are  parallel  to  the  external 
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field.  This  explains  why  computation  92  (L/K  =  -h)   is  even  slower  in 
convergence  than  computations  85  and  86  (l/k  =  -2).   It  is  significant  that 
the  existence  of  this  situation  is  strongly  brought  to  ones  attention  because 
of  the  rules  which  have  been  set  up  for  convergence  testing.  In  a  test  based 
on  the  examination  of  one  sequence  of  configurations  there  is  a  greater  chance 
that  one  would  fail  to  observe  this  near-critical  situation  since  the  sequence 
might  remain  entirely  in  one  set  of  configurations  during  the  computations. 
Furthermore^  it  should  be  noticed  that  the  standard  deviations  in  these  cases 
do  not  indicate  anything  unusual.  One  can  infer  from  the  small  standard 
deviations  that  once  the  convergence  conditions  were  satisfied  both  the  se- 
quences remained  in  the  one  class  of  states  which  was  most  probable. 

For  regions  in  which  series  expansions  of  the  long-range  order  and 
short-order  can  be  used  a  comparison  with  results  obtained  from  the  present 
Monte  Carlo  method  is  possible.  The  series  given  in  reference  (3)     have 
been  evaluated  for  a  few  cases  and  a  comparison  with  the  Monte  Carlo  results 
is  shown  in  Table  III. 

V.  THE  B0DY-CE1WEKED  CUBIC  LATTICE 

The  computations  for  this  lattice  are  not  as  extensive  as  those  for 
the  simple  cubic  lattice.  The  convergence  test  is  the  same  as  the  one  used 
with  the  simple  cubic  lattice  except  that  since  only  S.  and  f.(l)  are  calcu- 
lated during  each  iteration,  just  the  first  inequality  (8a)  is  required 
to  be  satisfied.  The  value  of  £..  was  set  at  0.02.   In  the  third  stage  of  the 
generation  of  the  configuration  sequence  (i«e.  after  the  convergence  condition 
was  satisfied)  the  results  for  the  two  independent  sequences  were  not  combined 
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but  instead  they  were  left   separate 0     Thus,   two  separate  estimates  of  the 
average  were  computed,   one  based  on  the  HP  sequence  and  the  other  based  on 
the  ET  sequence „     The  IIP  sequence  for  positive  K  was   started  with  an  initial 
configuration  in  which  all  of  the  spins  were  parallel  to  one  another  and 
parallel  to  the  external  field„     When  K  was  negative  the  initial  configuration 
for  the  IT  sequence  was  one  in  which  all  near est- neighbors  were  antiparallel. 
In  the  simple  cubic   lattice  computations  the  IT  sequence  always   started  with 
a  configuration  in  which  all  spins  were  parallel  to  each  other  and  to  the 
external  field,   even  when  K  was  negative . 

In  most  of  the  calculations   a  model  consisting  of  8  unit   cells 
on  an  edge,   and  therefore  1,024     sites,  was  used„     In  a  few  cases  a  smaller 
lattice  having  k  unit  cells  on  an  edge  was  used.     Twenty  seconds  was   required 
to  complete  one  iteration  on  the  larger  lattice „     This  program  was  not  quite 
as  efficient  in  its  use  of  computer  time  as  the  simple  cubic   lattice  program. 

The  sites  were  selected  for  detailed  consideration  in  a  systematic 
fashion.     The  two  sublattices  of  the   system  were  processed  separately;     that 
is,   all  the   "center"   sites  were  first  processed  sequentially,   then  all  "corner" 
sites  were  processed  sequentially  to  complete  one  iteration  for  one  lattice 0 

The  results  are  compiled  in  Table  IV  and  displayed  graphically  in 
Figs,   k  and  5°     In  Table  TV  the  upper  value  for  S  and  for  f(l)   is  obtained 
from  the  IT  sequence  and  the  lower  value  for  S  and  for  f (1)    is  obtained  from 
the  HT  sequence „     In  the  initial  calculations  the  larger  samples  with  M.  =  ^0 
and  ^M  -  50  were  computed,  but  to  conserve  on  computer  time  this  was   later 
reduced  to  M_  =  20  and^M  -  30„      In  all  cases  the  convergence  parameter  had 
the  value,    £     =  0»02o 
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Comparison  of  the  results   for  the   smaller  lattice,   having  four 
unit  cells  on  an  edge,  with  those  for  the  larger  lattice  shows  again  the 
relatively  strong  dependence  of  the  estimate  of  long-range  order  on  lattice 
size0     On  the  other  hand  the  short-range  order  estimate  again  is  not  so 
sensitive  to  change  in  lattice  size.     Comparison  of  results  for  the  in? 
sequences  with  those  for  the  HT  sequence  indicates  that  although  the  dif- 
ferences between  the  pairs  of  results  are  not   large,   they  are  frequently 
larger  than  the  spread  indicated  "by  the   standard  deviations „     This   is  not 
surprising,   since  the  correlation  that  exists  between  the  sequential  config- 
urations tends  to  reduce  the  standard  deviation  from  what  it  would  be   if  the 
configurations  in  the  chain  were  completely  independent „ 

VI.      CONCLUSION 

The  Monte  Carlo  method  appears  to  be  well  suited  to  the  computation 
of  short-range  order  in  an  Ising  lattice.      In  the  two-dimensional  lattice, 
where  the  results   can  be   checked  against  the  exact  treatment,   the  accuracy  of 
the  Monte  Carlo  results  for  a  20  x  20  array  is  quite  good.     With  one  exception, 
which  occurs   in  the  immediate  neighborhood  of  the  critical  temperature,    the 
errors  are    5$   or  less.     In  the  three  dimensional  systems,   arrays  of  ^12.  sites 
for  the  simple   cubic  and  1,02*1-  sites  for  the  body-centered  cubic  provided 
estimates  of  the  short-range  order  which  were  relatively  insensitive  to 
changes   in  the  size  of  the  array  even  in  the  neighborhood  of  the  critical 
temperature.     The  long-range  order,  on  the  other  hand,  was  found  to  be 
rather  sensitive  to  changes   in  size  of  the  array  near  the  critical  temperature. 
The  long-range  order  results   in  this   region  are  therefore  rather  crude.     An 
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accurate  estimate  of  long-range  order  near  the  critical  temperature  by  the 
present  method  does  not  seem  to  he  feasible  with  present  computing  equipment, 
"because  the  large  arrays  that  seem  to  he  required  demand  unreasonably  long 

computation  times. 

(3) 

As  pointed  out  by  Newell  and  Montroll,    the  primary  reason  for 

the  continued  interest  in  the  Ising  model  is  that  it  provides  a  simple  testing 
ground  for  theories  of  cooperative  phenomena.  The  present  results,  excepting 
perhaps  the  long-range  order  in  the  zero  field  case  near  T  ,  are  felt  to  he 
the  most  accurate  ones  now  available  and  they  have  been  tabulated  here  in 
considerable  detail  so  that  they  may  be  readily  compared  with  the  results  of 
other  methods.  Computations  using  a  simple  cubic  lattice  with  second  and 
third  neighbor  interactions  are  now  in  progress. 
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Figure  Captions 

FIGURE  1:   Long-range  order  S  in  the  simple  cubic  lattice 
shown  as  a  function  of  K  for  different  L  (solid 
curves)  and  L/K  (dashed  curves). 

FIGURE  2:   The  first  neighbor  short-range  order  parameter 
f (l)  in  the  simple  cubic  lattice  shown  as  a 
function  of  K  for  different  L  (solid  curves) 
and  L/K  (dashed  curves). 

FIGURE  3:   The  second  neighbor  short-range  order  parameter 
f (2)  in  the  simple  cubic  lattice  shown  as  a 
function  of  K  for  different  L  (solid  curves) 
and  L/K  (dashed  curves). 

FIGURE  k:        Long-range  order  S  in  the  body-centered  cubic 
lattice  shown  as  a  function  of  K  for  different 
L  (solid  curves)  and  L/K  (dashed  curves). 

FIGURE  5^  The  first  neighbor  short-range  order  parameter 

f(l)  in  the  body- centered  cubic  lattice  shown  as 
a  function  of  K  for  different  L  (solid  curves) 
and  L/K  (dashed  curves). 
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f(l 
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A 

2; 
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_0 
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Total 

0 

0.9919 

+ 

o.ooo4 

0.0079  + 

0.0003 

0.0080 

4 

0.0003 

50 

50 

156 

0 

0.9823 

+ 

0.0009 

0.0169  + 

0.0007 

0.0174 

4- 

0.0008 

25 

25 

55 

0 

0.9699 

+ 

0.0014 

0.0280  + 

0.0010 

0.0293 

4 

0.0011 

25 

25 

62 

0 

0.9424 

+ 

0.0011 

0.0516  + 

0.0006 

0.0549 

- 

0.0007 

100 

300 

403 

0 

0.8653 

4 

0.0060 

0.1073  + 

0.0021 

0.1181 

+ 

0.0026 

50 

50 

102 

0 

0.7^40 

+ 

0.0025 

0.1853  + 

0.0011 

0.2096 

+ 

0.0013 

100 

300 

4oo 

0 

0.6551 

+ 

0.0060 

0.2286  + 

0.0024 

0.2623 

+ 

0.0029 

50 

50 

161 

0 

0.6486 

4 

0.0024 

0.2313  + 

0.0010 

0.2652 

+ 

0.0013 

50 

50 

100 

0 

0.4976 

+ 

0.0090 

0.2841  + 

0.0024 

0.3315 

4 

0.0031 

50 

50 

116 

0 

0.3054 

+ 

0.0123 

0.3301  + 

0.0026 

0.3891 

4- 

0.0033 

50 

50 

100 

0 

0.1883 

4 

0.0091 

0.3405  + 

0.0016 

0.4013 

->- 

0.0020 

50 

50 

100 

0 

0.1603. 

+ 

o.oo¥) 

0.3713  + 

0.0007 

0.4344 

+ 

0.0008 

100 

300 

400 

0 

0.0781 

+ 

0.0062 

0.3708  + 

0.0013 

0.4341 

4 

0.0017 

50 

50 

100 

0 

0.0819 

+ 

0.0022 

0.4047  + 

0.0005 

0.4639 

+ 

0.0005 

100 

300 

4oo 

0 

0.0686 

4 

0.0018 

0.4211  + 

0.0004 

0.4755 

+ 

0.0003 

100 

300 

4oo 

0 

0.0563 

4 

0.0030 

0.4329  + 

0.0007 

0.481.8 

+ 

0.0006 

50 

50 

100 

0.0625 

0.8184 

+ 

0.0032 

0„l4it9  + 

0.0017 

0.1599 

+ 

0.0020 

50 

50 

102 

0.0625 

0.735^ 

+ 

o.oo44 

0.1958  + 

0.0021 

0.2192 

4- 

0.0025 

50 

50 

101 

0.0625 

0.5529 

+ 

0.0062 

0.2847  + 

0.0024 

0.3237 

* 

0.0028 

50 

50 

100 

0.0625 

0.3284 

+ 

0.0071 

0.3678  + 

0.0018 

0.4179 

4- 

0.0021 

50 

50 

1.00 

0.0625 

0.1834 

4 

0.0052 

0.4203  + 

0.0009 

0.4620 

+ 

0.0008 

50 

50 

100 

0.0625 

0.1025 

+ 

0.0086 

0.4709  + 

0.0009 

0.4946 

+ 

0.0007 

25 

25 

50 

0.0625 

O.0680 

+ 

0.0112 

0.5525  + 

0.0008 

0.4893 

+ 

0.0010 

25 

25 

50 

0.1250 

O.9601 

+ 

0.0010 

0.0372  + 

0.0007 

O.0388 

+ 

0.0008 

50 

50 

128 

0ol250 

O.9158 

+ 

0.0018 

0.0746  + 

Table   II 

0.0012 

O.0797 

4 

0.0014 

50 

50 

121 

-  29  - 


M 


K 

L 

s 

f(i 

i 

f( 

) 

_o 

4M 

Total 

0.2500 

0.1250 

0.8683 

+  0.0023 

0.1116  + 

0.0014 

0.1210 

+ 

0.0016 

50 

50 

105 

0.2381 

0.1250 

0.8383 

+  0.0022 

0.1340  + 

0.0013 

0.1459 

+ 

0.0015 

50 

50 

105 

0.2273 

0.1250 

0.8094 

+  0.0025 

0,1533  + 

0.0014 

0.1681 

+ 

0.0017 

50 

50 

105 

0.2174 

0.1250 

O.7678 

+  0.0054 

0.1788  + 

0.0024 

0.1972 

+ 

0.0028 

50 

50 

100 

0.2000 

0.1250 

0.6917 

+  0.0053 

0.2249  + 

0.0023 

0.2502 

+ 

0.0028 

50 

50 

100 

O.I667 

0.1250 

0.5303 

+  0.0050 

0.3084  + 

0.0019 

0.3445 

+ 

0.0022 

50 

50 

100 

0.1250 

0.1250 

0.3^15 

+  0.0045 

O.3880  + 

0.0013 

0.4286 

- 

0.0014 

50 

50 

100 

0.1250 

0.1250 

0.3369 

+  0.0019 

0.3892  + 

0.0006 

0.4297 

+ 

0.0007 

50 

50 

100 

0.0750 

0.1250 

0.21^5 

+  0.0060 

0.4424  + 

0.0010 

0.4722 

+ 

0.0010 

25 

25 

50 

0.1000 

0.1250 

0.0920 

+  0.0096 

0.5505  + 

0.0009 

0.4871 

+ 

0.0011 

25 

25 

50 

0.2000 

0.1250 

0.0613 

+  0.0023 

0.6266  + 

0.0023 

0.4342 

+ 

0.0031 

25 

25 

50 

0.2857 

0.2500 

0.9390 

+  0.0013 

0.0562  + 

0.0009 

0.0587 

r 

0.0010 

50 

50 

109 

0.2381 

0.2500 

0.8966 

+  0.0018 

O.0909  + 

0.0012 

0.0969 

+ 

0.0013 

50 

50 

103 

0.2000 

0.2500 

0.8127 

+  0.0024 

0.1548  + 

0.0014 

0.1670 

- 

0.0016 

50 

50 

102 

0.1250 

0.2500 

0.5655 

+  0.0034 

0.3082  + 

0.0015 

0.3341 

+ 

0.0016 

50 

50 

100 

0.0750 

0.2500 

0.3965 

+  0.0053 

0.3938  + 

0.0016 

0.4175 

+ 

0.0016 

25 

25 

50 

0.0625 

0.2500 

0.3678 

+  0.0050 

0.4091  + 

0.0013 

0.4303 

+ 

0.0013 

25 

25 

50 

0.0250 

0.2500 

0.2896 

+  0.0054 

0.4486  + 

0.0012 

0.4586 

-f 

0.0011 

25 

25 

50 

0.0500 

0.2500 

0.1964 

+  0.0038 

O.5056  + 

0.0007 

O.4798 

i 

0.0006 

50 

50 

100 

0.1000 

0.2500 

0.1566 

+  0.0038 

0.5394  + 

0.0007 

0.4785 

+ 

0.0007 

50 

50 

100 

0.1500 

0.2500 

0.1264 

+  o.oo4o 

0.5729  + 

0.0008 

0.4674 

+ 

0.0008 

50 

50 

100 

0.2000 

0.2500 

0.1030 

+  0.0037 

0.6171  + 

0.0013 

0.4366 

+ 

0.0015 

50 

50 

100 

0.2250 

0.2500 

0.0829 

+  0.0025 

O.6673  + 

0.0034 

0.3840 

+ 

o.oo4i 

25 

25 

61 

0.2500 

0.2500 

0.0564 

+  o.ooi4 

O.7960  + 

0.0024 

0.2274 

t 

0.0029 

50 

50 

122 

0.2750 

0.2500 

0.0423 

+  0.0019 

O.8655  + 

0.0023 

o.i46o 

4 

0.0026 

25 

25 

67 

Ta.ble  II (cont'd) 

30 


K 


f(l) 


f(2) 


M, 


0     AM       Total 


-0.3000 

-0.3500 

-0„5000 

0.1875 

0.1000 

0.0500 

-0.0250 

-0.1000 

-0.1500 

-0.1875 

-0.2500 

-0.3000 

-0.3750 

0.338 

0.2000 

0.0500 

-0.0500 

-OolOOO 
-0.1500 
-0.2000 
-0.2500 
-0.2750 
-0.3000 
-0.3250 
-0.3500 


0.2500 

0.2500 

0.2500 

0.7500 

0.7500 

0.7500 

0.7500 

0.7500 

0.7500 

0.7500 

0.7500 

0.7500 

0.7500 

1.250 

1.250 

1.250 

1»250 

1.250 

1.250 

1.250 

1.250 

1.250 

1.250 

1.250 

1.250 


0.0297 
0.0173 

o.oo64 

0.9473 

0.8424 
0.7^50 

0.5793 
0.4350 
0.3624 
0.3143 
0.2199 

0.1157 

0.0469 
0.9921 
0.9834 

0.9053 
0.7670 
0.6778 
0.5852 
0.5016 

0.4221 

0.3786 

0.3131 
0.2282 
0.1690 


+  0.0017 
+  0.0007 
+  0.0006 
+  o.ooi4 
+  0.0026 
+  0.0030 
+  o.oo4i 
+  0.0058 
+  o.oo64 
+  0.0059 
+  0.0029 
+  0.0021 

+  0.0014 

+  o.ooo4 
+  0.0007 
+  o.ooi4 
+  0.0020 
+  0.0023 
+  0.0027 
+  0.0033 
+  0.0035 
+  0.0061 
+  0.0022 

+  0.0033 
+  0.0020 


0.9033  + 
0.9504  + 
0.9876  + 
0.0498  + 
0.1405  + 
O.2165  + 
0.3381  + 

0.4377  + 
O.4915  + 

0.5331  + 
O.6561  + 

0.837^  + 
0.9408  + 
0.0078  + 
0.0163  + 
0.0902  + 
0.2103  + 
0.2840  + 

0.3578  + 
0.4261  + 
0.4960  + 
0.5384  + 
O.6216  + 
0.7358  + 
0.8088  + 


0.0020 
0.0010 
0.0006 
0.0010 

0.0016 
0.0019 

0.0015 

0.0017 
0.0018 

0.0021 

0.0037 
0.0028 

0.0012 

0.0003 
0.0005 
0.0019 

0.0012 
0.0011 
0.0011 
0.0012 
0.0013 
0.0027 
0.0026 
0.0032 
0.0020 


0.1040  + 
0.0519  + 
0.0126  + 
0.0513  + 
0.1452  + 
0.2216  + 
0.3324  + 

o.4ooi  + 
0.4206  + 

0.4223  + 

0.3416  + 
0.1613  1 
0.0589  + 
0.0079  + 
0.0165  + 
0.0911  + 
0.2065  + 
0.2697  + 
0.3239  + 
0.3598  + 
0.3754  + 
0.3679  + 
0.3169  + 
0.2280  + 
0.1711  + 


0.0023 

0.0011 

0.0006 

0.001.1 

0.0017 
0.0019 
0.0016 
0.0013 
0.0015 
0.0016 
0.0038 
0.0030 

0.0012 
0.0003 
0.0005 
0.0014 
0.0012 
0.0011 
0.0009 
0.0008 
0.0008 
0.0021 
0.0023 
0.0028 
0.0018 


25 
50 
25 
25 
25 
25 
25 
25 
25 
25 
25 
25 
25 
50 
50 
50 
50 
50 
50 
50 
50 
25 
50 
25 
50 


25 
50 
25 
25 
25 
25 
25 
25 
25 
25 
25 

25 

25 
50 
50 
50 
50 
50 
50 
50 
50 
25 
50 
25 
50 


81 
106 

79 

62 

52 

51 

50 

50 

50 

50 

55 

89 

67 

101 

101 

100 

100 

100 

100 

100 

100 

50 

129 

52 

116 


M.. 


K 

L 

s 

f( 

1J 

1 

f( 

I. 

[ 

_0_ 

AM 

Total 

■0.4000 

1.250 

O.O96I 

+ 

0.0020 

0.8940 

4 

0.0017 

0.0993 

+ 

0.0015 

25 

25 

58 

0.5000 

1.250 

0.0298 

+ 

o„ooo7 

0.9686 

+ 

0.0006 

0.0307 

+ 

0.0006 

50 

50 

134 

0.2000 

0.4000 

0.8834 

+ 

0.0022 

0.1035 

4 

0o0015 

0.1091 

+ 

0.0017 

25 

25 

64 

0.1750 

0.3500 

0.8171 

+ 

0.0033 

0.1543 

4 

0.0018 

0.1645 

+ 

0.0021 

25 

25 

79 

0.1500 

0.3000 

0.7121 

+ 

0.0047 

0.2251 

4 

0.0025 

0.2430 

+ 

0.0028 

25 

25 

95 

0.0500 

0.1000 

0.1080 

+ 

0.0123 

0.5233 

+ 

0.0008 

0.4964 

+ 

0.0009 

25 

25 

50 

0.1000 

0.2000 

0.1364 

4 

0.0086 

0.5443 

+ 

0.0011 

0.4821 

t 

0.0012 

25 

25 

50 

0.1500 

0.3000 

0.1544 

+ 

0.0076 

0.5687 

4 

0.0013 

0.4675 

4 

o.ooi4 

25 

25 

50 

0.2000 

0.4000 

0.1604 

+ 

0.0052 

0.6035 

+ 

0.0018 

0.4341 

+ 

0.0015 

25 

25 

51 

0.2500 

0.5000 

0.1260 

+ 

0.0025 

0.7397 

+ 

o.oo4o 

0.2806 

4 

0.0049 

25 

25 

60 

0.3000 

0.6000 

0.0863 

+ 

0.0017 

0.8622 

4 

0.0021 

o.i.4o4 

+ 

0.0023 

25 

25 

87 

0.0500 

0.2000 

0.1679 

+ 

0.0084 

0.5133 

+ 

0.0012 

0.4871 

4- 

0.0011 

25 

25 

50 

0.1000 

0.4000 

0.2468 

+ 

O.OO67 

0.5182 

t 

0.0016 

0.4629 

+ 

0.0013 

25 

25 

50 

0.1500 

0.6000 

0.2915 

+ 

0.0059 

0.5245 

+ 

0.0018 

0.4410 

4 

0.0013 

25 

25 

50 

0.2000 

0.8000 

0.3204 

+ 

0.0055 

0.5362 

+ 

0.0020 

0.41.68 

4 

0.0015 

25 

25 

50 

0.2500 

1.0000 

0.3260 

+ 

0.0055 

0.5665 

4 

0.0024 

0.3823 

+ 

0.0017 

25 

25 

50 

0.3000 

1.2000 

0.2889 

+ 

0.0030 

0.6466 

4 

0.0039 

0.3046 

+ 

0.0036 

25 

25 

ill 

0.0500 

0.4000 

0.3023 

+ 

0.0052 

0.4759 

+ 

o.ooi4 

0.4536 

4 

0.0013 

25 

25 

50 

0.1000 

0.8000 

0.4638 

+ 

0.0051 

0.4226 

+ 

0.0015 

0.3886 

+ 

0.0012 

25 

25 

50 

0.1500 

1.2000 

0.5635 

+ 

0.0051 

0.3736 

+ 

0.0018 

0.3357 

+ 

0.0014 

25 

25 

50 

0.2000 

1.6000 

0.6321 

+ 

0.0042 

0.3313 

+ 

0.0016 

0.2958 

4 

0.0012 

25 

25 

50 

0.2500 

2.0000 

0.6898 

+ 

0.0038 

0.2906 

4- 

0.0015 

0.257^ 

4 

0.0014 

25 

25 

50 

0.3000 

2.4000 

0.7282 

4 

0.0033 

0.2595 

4- 

0.0014 

0.2314 

4- 

0.0012 

25 

25 

50 

Table  II (cont'd) 
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K  (L  =  0)  S  f(l) 


Monte 

Monte 

Carlo 

Series 

Carlo 

Series 

0.5000 

0.9919 

0.99^5 

0.0079 

0.0054 

0.4000 

0.9699 

0.9795 

0.0280 

0.0195 

0.3380 

0.9424 

0.9504 

0.0516 

0.0453 

0.2000 

— 

— 

0.3713 

0.3756 

O.I665 

— 

— 

0.4048 

0.4050 

0.1430 

-- 

— 

0.4210 

0.4217 

Table   III 
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K 


f(D 


M, 


AM         Total 


0.3 


0.225 


0.225 


0.1875 


0.175 


0.175 


0.1625 


0.15 


0.1375 


0 


0 


0.9804  +  0.0007 
0.9800  +  0.0018 

0.8966  +  0.0060 
0.8931  +  0.0129 

0.9067  +  0.0019 
0.9032  +  0.0032 

0.8023  +  0.0028 
0.7850  +  0.0095 

0.7127  +  0.0109 

O.6717  +  0.0202 

0.6884  +  0.0050 
O.6697  +  0.0102 

0.5216  +  0.0091 
0.5246  +  0.0103 

0.2159  +  0.0132 
0.1602  +  0.0101 

0.0938  +  0.0095 
0.1091  +  0.0071 


O.OI89  +  0.0007 
0.0188  +  0.0012 

0.0908  +  0.0045 
0.0866  +  o.oo64 

0.0831  +  0.0015 
0.0847  +  0.0016 

0.1607  +  0.0018 
0.1662  +  0.0034 

0.2173  +  0.0063 
0.2251  +  0.0088 

0.2251  +  0.0027 
0.2321  +  0.0035 

0.2981  +  0.0033 
0.2982  +  0.0032 

0.3734  +  0.0027 
0.3850  +  0.0017 

0.4030  +  0.0023 
0.4099  +  0.0013 


50 


50 


50 


50 


50 


50 


50 


50 


50 


50 


50 


50 


50 


50 


50 


50 


50 


101 


100 


103 


101 


100 


50    109 


100 


101 


100 


Table  IV 
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K 


0.125 


0.1 


s 

f(i) 

_2 

AM 

Total 

0.0781*  +  0.0071 
0.0684  +  0.0051 

0.4179  +  0.0017 

0.4223  +   0.0013 

50 

50 

100 

O.1366  +  0,0103 
O.1433  +  0.0107 

0.4430  +  0.0030 
0.4376  +  0.0029 

50 

50 

100 

0.1 


0.05 


0.05 


0.2 


0.1625 


0.15 


0.125 


0.1 


0 


0 


0.05 


0.05 


0.05 


0.05 


0.05 


0,o44i  +  0.0038 
0.0484  +  0.0036 

0.1031  +  0.0098 
O.1169  +  O0O074 

0.0453  +  0.0062 
0.0318  +  0.0024 

O.8809  +  0.0029 
0.8823  +  0.0031 

O.7136  +  0.0056 
0.7137  +  0.0052 

O.6190  +  0.0093 
O.6018  +  0.0073 

0.3135  +  0.0102 
O.3165  +  0.0130 

O.1530  +  0.0076 
0.1516  +  0.0094 


0.4404  +  0.0012 
0.4430  +  0.0010 

0,4635  +  0.0036 
0.4714  +  0.0022 

0.4684  +  0.0026 
0.4767  +  0.0018 

0.1044  +  0.0023 
0ol030  +  0.0024 

0.2190  +  0.0032 
0.2185  +  0.0030 

0.2715  +  0.0042 

0.2769  +  0.0036 

0.3869  +  0.0027 
0.3853  +  0.0039 

0.4351  +  o.ooi4 
0.4330  +  0.0023 


50 


50 


50 


20 


20 


20 


20 


20 


50 


50 


50 


30 


30 


30 


30 


30 


100 


100 


101 


58 


56 


50 


50 


50 


Table  IV (cont'd) 
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K 

0.15 


0.125 


0.1 


0.05 


0.15 


0.125 


0.1 


0.05 


0.175 


0.15 


0.15 


0.15 


0.15 


0.15 


0.25 


0.25 


0.25 


0.25 


0.35 


0.35 


s 

f(i) 

_H 

AM 

Total 

0.7657  +  0.0047 
0.7663  +  0.0047 

0.1906  +  0.0031 
0.1902  +  0.0028 

20 

30 

51 

0.6051  +  0.0062 
0.6113  +  0.0063 

0.2857  +  0.0034 
0.2836  +  0.0030 

20 

30 

53 

0.4388  +  0.0077 
0.4441  +  0.0086 

0.2338  +  o.oo64 
0.2349  +  0.0066 

0.8371  +  0.0034 
0.8338  +  0.0034 

0.7361  +  0.0047 

0.7379  +  0.0042 

0.6084  +  0.0052 
0.6056  +  0.0053 

0.3789  +  0.0055 
0.3712  +  0.0069 

0.9222  +  0.0037 
O.9198  +  0.0020 

0.8718  +  0.0033 
0.8795  +  0.0027 


0.3652  +  0.0031 
0.3632  +  0.0033 

0.4517  +  0.0018 
0.4472  +  0.0016 

0.l4l6  +  0.0025 
0.1430  +  0.0024 

0.2130  +  0.0032 
0.2119  +  0.0027 

0.2905  +  0.0029 
0.2946  +  0.0030 

0.4089  +  0.0023 

0.4108  +  0.0024 

0.0720  +  0.0030 
0.0741  +  0.0018 

0.1143  +  0.0026 
0.1081  +  0.0023 


20 


20 


20 


20 


20 


20 


20 


20 


30 


30 


30 


30 


30 


30 


30 


30 


50 


50 


51 


50 


50 


50 


50 


50 


Table  IV( cont'd) 
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K 


0.125 


0.1 


0.05 


-0.05 


-0.1 


-0.15 


-0.2 


-0.25 


-0.05 


-0.15 


0.35 


0.35 


0.35 


0.35 


0.35 


0.35 


0.35 


0.35 


0.5 


0.75 


0.8135  +  0.003^ 
O.8076  +  0.0040 

0.7105  +  0.0042 

0.7164  +  o.oo4o 

0.4970  +  0.0051 
0.4943  +  0.0049 

0.2451  +  0.0060 
0.2463  +  O.OO85 

0.1849  +  0.0050 
0.1870  +  0.0059 

0.1387  +  0.0033 
0.1360  +  0.0028 

0.0581  +  0.0024 

0.0646  +  0.0021 

0.0296  +  0.0011 

0.0268  +  0.0015 

0.3431  +  0.0057 
0.3467  +  0.0065 

0.3060  +  0.0054 
0.3044  +  0.0036 


f(l) 

0.1596  +  0.0026 
0.1649  +  0.0029 

0.2329  +  0.0027 
0.2293  +  0.0026 

0.3621  +  0.0025 

0.3613  +  0.0024 

0.4922  +  0.0023 
0.4916  +  0.0018 

0.5304  +  0.0020 

0.5361  +  0.0017 

0.5954  +  0.0021 

0.6010  +  0.0034 

0.8300  +  0.0050 
0.8501  +  0.0028 

0.9386  +  0.0047 
0.9461  +  0.0017 

0.4591  +  0.0025 
0.4601  +  0.0157 

0.5278  +  0.0032 
0.5321  +  0.0167 


M 


20 


20 


20 


20 


20 


20 


20 


20 


20 


20 


AM       Total 


30  50 


30  50 


30  50 


30  50 


30  50 


30  51 


30  63 


30  51 


30  50 


30  50 


Table  IV (cont'd) 
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K 


0.05 


-0.1 


-0.15 


-0.2 


-0.25 


-0.3 


-0.25 


-0.2 


1.0 


1.0 


1.0 


1.0 


1.0 


1.0 


1.25 


2.0 


S 

0.8729  +  0.0025 
0.8725  +  0.0025 

0.5127  +  0.0067 
0.5133  +  0.0088 

o.it-107  +  0.0056 
0.4119  +  0.0074 

0.2704  +  0.0035 
0.2734  +  0.0028 

0.1299  +  0.0038 
0.1222  +  0.0025 

0.0655  +  0.0066 
0.0548  +  0.0016 

0.2003  +  O.OO78 
0.1921  +  0.0025 

0.6734  +  0.0057 
0.6728  +  0.006k 


f(l) 

12 

Am 

Total 

0.1174  +  0.0021 
0.1180  +  0.0022 

20 

30 

55 

0.39^7  +  o.oo44 
0.3945  +  0.0063 

20 

30 

50 

0.4715  +  0.0036 
0.4738  +  o.oo44 

0.64-72  +  0.0055 
0.6472  +  0.0037 

0.8472  +  0.0054 
0.8604  +  0.0028 

0.9247  +  0.0091 
O.9409  +  0.0016 

O.7744  +  0.0091 
O.7877  +  0.0028 

O.2991  +  0.0051 
0.3002  +  0.0055 


20 


20 


20 


20 


20 


20 


30 


30 


30 


30 


30 


50 


56 


52 


50 


50 


30     50 


Table  IV (cont'd) 
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1  (L  =  0)  S  f(l) 


Monte 

Carlo  Series 


0.300  0.9804  O.9806 

O.98OO 

0.150 


0.050 


Table  V 


Monte 
Carlo 

Series 

O.OI89 

0.0188 

0.0188 

0.373^ 

0.3955 

0.3850 

0.4684 

0.4742 

O.4767 
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